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Abstract 


In  this  report rw?  investigateJ^he  smoothing  problem  for  two-point 
boundary  value  models.  We  show  thatT by  triangularizing  the  relevant 
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Hamiltonian  system,  ws  can  derive  a  new  smoothing  algorithm  which  is 

twice  as  fast  as  the  one  which  has  been  obtained  by  diagonalizing  the 
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1.  Introduction 


When  modeling  physical  systems  with  a  spatial  independent  variable, 
the  a  priori  information  is  usually  given  in  the  form  of  boundary  conditions 
rather  than  initial  conditions.  Consider  the  following  linear,  two-point 
boundary  value  model  for  to  <  t  <  tj : 

x(t)  =  Ax[t) + Bu[t)  (1) 

r  =  V0x[t0)  +  V1x{t1)  (2) 

where  u  is  a  zero-mean,  unit  intensity  white  noise  with  m  components,  r  is  a 
zero-mean  random  vector  with  n  components  and  covariance  Q,  and  is 
uncorrelated  with  u,  x  is  the  n-component  state  vector,  and  A,  B,  Vq,  V1  are 
appropriate  matrices.  We  have  continuous  measurements  on  the  interval 
[to,  tj  given  by 

y(t)  =  Cx(t)+u(t) 

where  v  is  a  zero-mean,  unit  intensity  white  noise  with  p  components  and  is 
uncorrelated  with  u  and  r. 

As  is  customary,  we  will  assume  that  the  model  is  well-posed;  i.e.,  that 
u  and  r  give  rise  to  a  unique  x.  This  will  be  the  case  if  and  only  if  the  matrix 
V0  +  is  nonsingular,  in  which  case  we  can  assume,  without  loss  of 

generality,  that 

V0  +  V1<P(tl,t0)  =  I 

where  <P  is  the  state  transition  matrix  of  (1). 

The  problem  of  developing  stable  recursive  algorithms  for  the  linear 
least-squares  smoothed  estimate  of  the  state  x  given  the  measurements  y. 
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and  for  the  associated  error  covariance,  has  been  addressed  by  several 
authors  [l]-{5].  In  these  references,  the  Hamiltonian  equations 
characterizing  the  smoothed  estimate  are  solved  by  the  use  of  a 
diagonalizing  Riccati  transformation.  This  transformation  requires  the 
solution  of  two  Riccati  equations.  One  can  also  solve  the  Hamiltonian 
equations  by  triangularizing  them,  a  procedure  that  requires  the  solution  of 
only  one  Riccati  equation.  Although  triangularization  was  mentioned  in  [1] 
as  an  option,  the  necessary  equations  and  a  comparison  with  diagonalization 
were  not  provided.  It  turns  out  that  triangularization  leads  to  substantial 
savings  in  computation.  In  fact,  as  we  will  show,  the  number  of  required 
floating  point  operations  is  cut  in  half. 

2.  Triangularization  and  Shooting 

The  Hamiltonian  equations  for  the  smoothed  estimate  x  are  derived  in 
various  equivalent  forms  In  [l]-[5].  The  form  we  will  use  Is  given  below: 


Now  let  N  be  the  solution  to  the  following  Riccati  equation: 

N[t)  =  -N[t)A - A'N[t)  +  N[t)BB'N[t) - C'C,  N(t,)  =  0  (5) 

The  coordinate  change 

\m]j  i  ow 

/law. 
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transforms  eqns.  (3) -(4)  into 


x(£)  _  A-BB'j 

o 


'N(t)  BB' 
-( A-BB 


'  m  r  o  ] 

mt)Ylp[t)\+[-Cfylt)\  t€  0,1 


v0+9N(t0)  -9p(t0)i  rvj  Qo%,t0)  p^)]  m 

_  -V{N{t0)  V/iptto)]  L°  f-Vi^.tojJLpUJJ 

Eqns.  (6)-(7)  constitute  a  triangularized  two-point  boundary  value  problem 

which  can  be  solved  by  a  standard  shooting  procedure.  Of  course,  shooting 

could  be  applied  directly  to  eqns.  (3) -(4),  but  the  resulting  system  would  not 

be  stable.  The  point  of  using  a  Riccati  transformation  is  to  produce  stable 

differential  equations.  If  [A,  B)  is  stabilizable  and  (A,  C)  is  detectable,  then 

(A-BB'N)  will  be  stable  in  the  direction  of  increasing  t. 

To  solve  eqns.  (6) -(7),  first  solve  the  following  equations: 

*0(£)  =  (A ~ BB'N{t))x0(t)  +  BB'p0,  x0[t0)  =  0  (8 ) 

p0 (t)  =  -(A - BB'N(t))'p0 it) - C'y(t),  p0(tj)  =  0  ( 9 ) 

Then  it  can  easily  be  checked  by  differentiation  that 


m = x0it)  +  nt,t0m0)+ mm  y'mmpM) 

p(t)=p0(t)+  ¥''(£j,t)p(ti) 

where  ‘Pis  the  state  transition  matrix  of  ( A-BB'N)  and 


(ID 


M(t)  =  {A - BB'N{t))M[t )  +  M{t)[A - BB'N[t)Y  +  BB',  M(t0)  =  0  ( 1 2) 

Plugging  xitj  and  p(^)  obtained  from  eqns.  (10)-(11)  into  eqn.  (7)  gives 

r-9  vq  - mt0 )  ^.toj-^.tojipctonpPoCtoj-vixo^n 

[v{  olntM-GitM  mm)  JjLpftjJ  l  -v^o(to)  J 

The  solution  of  this  equation  is  then  used  in  eqn.  (10)  to  yield  the  smoothed 
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state  estimate.  The  well-posedness  of  eqns.  (l)-(2)  and  the  uniqueness  of 
linear  least-squares  estimates  guarantee  that  the  coefficient  matrix  in  eqn. 
(10)  is  nonsingular. 


3.  Error  Covariance 

The  smoothing  error  x  satisfies  [1],  [51: 


r*(t) 

"  A 

BB'Tx(t)  ‘ 

i 

’Bu(t)" 

1 

1 

. J 

C'C 

-A'JL-A(t) 

T 

_C'u(£) 

.  telto.tj 


X  -Sptto)  ' 

1 

Vj  Q0%,to)  Txft)  ‘ 

Jr] 

.0  VJ'l-Att,). 

~r 

L°J 

(14) 


(15) 


After  the  coordinate  transformation 


x(t)‘ 

1 _ 

°1 

x(t)  ' 

.set). 

Nit) 

1 

-A(t). 

eqns.  (14)-(15)  become  triangularized  as  follows: 


x(t) 

~A-BB'N[t) 

BB' 

’m 

x 

Bu(t) 

.let). 

0 

-{A-BB'N{t))\ 

.m. 

*r 

N{t)Bu{t)  +  C'v{t)\ 

' v0+QN[t0 )  -g- 

’xltj 

i 

V,  Q4>%,t0)  ' 

*(tx)' 

V 

~V{N(t0)  V{ 

T 

0  J-V^'Mo). 

0 

L  t  e  Mil  (16) 


(17) 


To  proceed,  we  again  use  shooting.  As  before,  we  can  write 

x(t)  =  x0(t)+‘P(t,to)x(t0)  +  M(t)‘P'(t1.t)^(t1) 

where 

x0(t)  =  (A  -  BB'JV(t))x0(t)  +  BB'$o (t)  +  Bu(f),  x0(tQ)  =  0 


(18) 

(19) 

(20) 
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$0(t)  =  ~{A  - BB'N(t)YS0{t)  +  N{t)Bu{t)  +  C'vit),  Z0lt1)  =  0  (21) 

Plugging  x(t,)  and  £(£o)  obtained  from  eqns.  (18)-(19)  into  eqn.  (17)  gives 


*{t0) 

uw. 


=N+wP°(t>>l 

W  LfoM 


where  F  is  the  coefficient  matrix  in  eqn.  (13)  and 

W  = 

As  a  result,  we  can  rewrite  eqn.  (18)  as 


-Vi  0 ' 
0  -V' 


x(t)  =  xJt)+R(t)F 


-i 


+  W 


x0[tj 

5o(«o) 


where 


(22) 


R{t)  =  [nt,t0) 

We  are  now  in  a  position  to  determine  the  error  covariance  from  eqns.  (20)- 

(22). 

First,  using  standard  techniques,  one  can  show  that  if 


£{#(t)  =  E[£0(t)£(t)] 


then 


tig  (t)  =  -(A  -  BB'N[t))'Z(o  (t)  -  Zfo  (t)(A  -  BB'N{t))  -  N[t)BB'N[t)  -  C'C,  S(o  (t, )  =  0 

But  in  fact,  since  their  difference  satisfies  a  homogeneous  linear  differential 
equation  with  zero  initial  condition, 

and  thus 
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(23) 


E[£0(t)£'(s)]  = 


f‘f"(s.t)N(s),  t  <  s 
[N{tmt,s),  t> s 


It  is  shown  in  [6]  that  the  process  u[t)+B'^0[t)  is  white.  Therefore,  if 

Zio(t)  =  E[x0{t)x'0{t)] 

then,  referring  to  eqn.  (12), 

and 


E(x0(t1)x'(t)i=  n^.twm 


One  can  also  check,  using  eqns.  (20)-(21).  (23)  that 

E[x0[t)Z'0{t0)]  =  0 

Taking  the  covariance  of  both  sides  of  eqn.  (22),  we  get 


P(t)  =  E[*(t)*'(t)] 

=  M(t)  +  P(t)F~1 


-R[t)F 


-I 


Q+QmtJQ+VMtN  -QNW 

-VWtJQ  V{N(t0)V i 

Vl  «F(ti ,  t)M(t)  _  * t)y,  ojF'-'P'It) 


F'-'R'it) 


(24) 


4.  Conclusions 

We  have  derived  a  new  smoothing  algorithm  for  two-point  boundary 
value  models  based  on  triangularizing  the  Hamiltonian  system.  If  one 
carefully  tabulates  the  matrix  multiplications  needed  to  determine  all 
intermediate  quantities,  one  finds  that  for  each  point  t  at  which  the 
smoothed  estimate  and  its  error  covariance  must  be  computed,  18n3  floating 
point  operations  are  required  (terms  of  order  n2  or  less  have  been  ignored). 
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I  ■ 1  1  -  -  ■  ' —  -  "  -  '  -  -  --  -  . 

Alternatively,  for  the  algorithm  detailed  in  [1],  38n3  floating  point  operations 
are  required.  Our  algorithm  is  therefore  slightly  more  than  twice  as  fast. 
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